home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Mac Magazin/MacEasy 19
/
Mac Magazin and MacEasy Magazine CD - Issue 19.iso
/
Wissenschaft & Technik
/
FFTs for RISC 1.1
/
fftTest.c
< prev
next >
Wrap
C/C++ Source or Header
|
1996-01-30
|
3KB
|
144 lines
/* A program to test and time complex forward and inverse fast fourier transform routines */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "fftlib.h"
#define MACTIMING 1 /* true means time fft with macintosh calls */
#if MACTIMING
#include <timer.h>
#endif
#define NUMROWS 1 /* process Matrix of NumRows different ffts of length N */
#define N 1024 /* size of FFT must be a power of 2 */
#define NTIMES 20 /* number of timings,invalid if too big (if a[0][0].Re = 0|nan)*/
typedef struct{
float Re;
float Im;
} Complex;
void main(){
float *Utbl;
Complex (*a)[N];
#if MACTIMING
UnsignedWide TheTime1;
UnsignedWide TheTime2;
UnsignedWide TheTime3;
double TheTime;
#endif
long i, il;
long TheErr;
long M;
Utbl = (float *) calloc((N/4+1),sizeof(float));
if (Utbl==0)
TheErr = 2;
else
TheErr = FFTInit(&M, N, Utbl);
if(!TheErr){
a = (Complex (*)[N]) calloc(NUMROWS*N,sizeof(Complex));
if (a == 0) TheErr = 2;
}
if(!TheErr){
/* set up a simple test case */
for (il=0; il<NUMROWS; il++){
for (i=0; i<N; i++){
a[il][i].Re = sqrt(il+i+.77777);
a[il][i].Im = (il+i+.22222)*(il+i+.22222) / N - N/2;
}
a[il][0].Re = N+3;
a[il][1].Re = 1-N;
}
#if MACTIMING
Microseconds(&TheTime1);
for (i=0;i<NTIMES;i++){ /* do NTIMES times for timing */
ffts((float *)a, M, NUMROWS, Utbl);
}
Microseconds(&TheTime2);
for (i=0;i<NTIMES;i++){ /* do NTIMES times for timing */
iffts((float *)a, M, NUMROWS, Utbl);
}
Microseconds(&TheTime3);
TheTime = (double)(TheTime2.hi - TheTime1.hi) * 65536.0 * 65536.0;
TheTime = (TheTime + (double)(TheTime2.lo - TheTime1.lo))/1000.0;
printf("Time fft = %12f ms.", TheTime/NTIMES/NUMROWS, a[0][0].Re);
TheTime = (double)(TheTime3.hi - TheTime2.hi) * 65536.0 * 65536.0;
TheTime = (TheTime + (double)(TheTime3.lo - TheTime2.lo))/1000.0;
printf(" ifft = %12f ms. a[0][0].Re= %6e\n", TheTime/NTIMES/NUMROWS, a[0][0].Re);
#else
printf("start timing \n");
for (i=0;i<NTIMES;i++){ /* do NTIMES times for timing */
ffts((float *)a, M, NUMROWS, Utbl);
iffts((float *)a, M, NUMROWS, Utbl);
}
printf("end timing \n");
#endif
printf("\n");
/* set up a simple test case */
for (il=0; il<NUMROWS; il++){
for (i=0; i<N; i++){
a[il][i].Re = sqrt(il+i+.77777);
a[il][i].Im = (il+i+.22222)*(il+i+.22222) / N - N/2;
}
a[il][0].Re = N+3;
a[il][1].Re = 1-N;
}
ffts((float *)a, M, NUMROWS, Utbl);
if (N*NUMROWS <= 256){
for (il=0; il<NUMROWS; il++){
printf("atrans = [ \n");
for (i=0; i<N; i++)
printf(" %+20.15e + j * ( %+20.15e ) \n", a[il][i].Re, a[il][i].Im);
printf("]; \n");
}
}
else { /* abbreviate big output */
printf("the first fft's last 32 values are: \n");
for (i=N-32; i<N; i++)
printf(" %+20.15e + j * ( %+20.15e ) \n", a[0][i].Re, a[0][i].Im);
}
iffts((float *)a, M, NUMROWS, Utbl);
if (N*NUMROWS <= 256){
for (il=0; il<NUMROWS; il++){
printf("\n aitrans = [ \n");
for (i=0; i<N; i++)
printf(" %+20.15e + j * ( %+20.15e ) \n", a[il][i].Re, a[il][i].Im);
printf("]; \n");
}
}
else { /* abbreviate big output */
printf("\n the first ifft's last 32 values are: \n");
for (i=N-32; i<N; i++)
printf(" %+20.15e + j * ( %+20.15e ) \n", a[0][i].Re, a[0][i].Im);
}
free (a);
free (Utbl);
return;
}
else
if(TheErr==2) printf(" out of memory ");
else printf(" error ");
return;
}